Beam propagation in finite size photonic crystals and metamaterials 
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The recent interest in the imaging possibilities of photonic crystals (superlensing, superprism, 
optical mirages etc..) call for a detailed analysis of beam propagation inside a finite periodic 
structure. In this paper, an answer to the question "where does the beam emerge?" is given. 
Contrarily to common knowledge, it is not always true that the shift of a beam is given by the 
normal to the dispersion curve. This phenomenon is explained in terms of evanescent waves and a 
£\j I renormalized diagram that gives the correct direction is given. 
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PACS numbers: 42.70 Qs, 42. 25. Fx 



INTRODUCTION AND SETTING OF THE PROBLEM 

Some beautiful experiments and numerical works have shown that it was possible to obtain quite unusual behaviors 
of light propagation inside meta-materials and photonic crystals (PhCs) [EH, S S H 0> @| • In particular, the near 
field properties of meta-materials are intensively studied in view of the possibility of designing superlenses |T(| or 



cloaking devices 



In these structures, the evanescent waves play an important role: the point of this work is to to quantify the 
importance of the evanescent waves and to give a theoretical insight into the propagation of a beam [12j | inside a 
finite-size PhC. In principle, the group velocity (13 , allows to determine where the beam emerges from the PhC 
O . (see figure 1) by computing the shift A. We show that, in finite-size structures, the shift is not always correctly 

predicted by the normal to the isofrequency curves. 
-H . This fact is due to two reasons. First, in finite-size structures, there are evanescent waves near the boundaries 
£> ' which can contribute to the propagation of the beam . Therefore the field inside the structure comprises not only 



Bloch modes but evanescent waves as well [22| . The latter can have a strong influence on the behavior of the beam. 
Second, due to multiple scattering, the emerging field is a sum of beams that can strongly interfere. If the beams 
. are well separated spatially, on can clearly distinguish where the first beam emerges from the structure. If the beams 
overlap strongly, the resulting field can be strongly deformed, and it becomes difficult to define an "emerging point". 
Throughout this work, we use time-harmonic fields, with a time-dependence of exp(— i0t). The fields are assumed 
qq ' to be z-independent (this is the direction of invariance of the photonic crystal). The vectorial diffraction problem 
can then be reduced to the study of the two usual cases of polarization: s-polarization (electric field parallel to z) 
or p-polarization (magnetic field parallel to z). The wavenumber in the medium surrounding the photonic crystal is 
•i-H , denoted by k . The incident field is a limited beam whose z component is given by 

03 ; u l (x,y) = [ ° A{k x )e i{k * x+k v ov) dk x (1) 

J-ko 



where: k y0 = y/k^ — k% and A(k x ) is the spectral amplitude (for instance, it can be chosen gaussian: A(k x ) = 

e - : V( fc ^~ fc "i) 2 ) ; where k rn = kosmO, and is the mean angle of incidence of the beam and w its waist. The point 
where the incident beam enters the photonic crystal is defined as the barycenter of the beam. Its abscissa is given 

f x |ia 1 (a:,0) J dx 



by: Xi = ...... ,,.ir 



2 



SHIFT OF THE FIRST TRANSMITED BEAM 



The crystal is described as a stack of gratings (figure 1) and we assume that in the spectral domain defined by the 
above beam, the ratio between the wavelength and the period d of the gratings is such that there is only one reflected 
and one transmitted order (the condition ko < tr/d is sufficient). Given this hypothesis, the reflected and transmitted 
fields can be expressed as: 



u r {x, y) = J A (k x ) r N (k x ) e^ x ~ k ^dk x (2) 
u t (x, y) = [ A (k x ) t N (k x ) e l ^ x+k ^dk x (3) 



where and are the reflection and transmission coefficients. For a given k Xl there exists a unique real 2x2 
matrix T^r 16|, 17 1 with determinant 1, satisfying the following relation: 



- N 



1 + rN A M,, 1 ) (4) 



ikyo (1 - tn) J \ ikyo 



It is the dressed transfer matrix of the total structure [la]. Let us denote by 7 and 7 _1 the eigenvalues of Tat 
and by v = (</>n, </>2i) , w = (0i2, ^22) the associated eigenvectors (Tv = 7V, Tatw = 7~ 1 w). The reflection and 
transmission coefficients can be written in the following form: 

(7 2 - 1) / . 7(1 -ff- 1 /) _ 

where, denoting q(x, y) — {ik y0 y — x) / '(ik y oy + x) , the functions / and g are defined by g (k, 9) = q (v) , / (k,6) = q (w) 
and v is chosen such that \g\ < 1 in the conduction bands. 
The following expansions hold: 

+00 

r N (k,6) = g + (g-f)J2^ P \9\ 2p (6) 
P =i 

t N (k,6) = (l-\g\ 2 hJ2^ P \9\ 2p (7) 

p=0 

These series expansion show that, due to the multiple scattering inside the photonic crystal, the transmitted (and 
reflected) fields consist of a sum of beams Ut(x, y) — Ylp=o u ti where: 

uf (x, y)= [ A(k x )(l - \g\ 2 ) \g\ 2p ^e^+^dk, (8) 



the first transmitted field: A = X t — X i: where: X t = n \ 



The position where a beam emerges from the photonic crystal is defined as its barycenter. We are interested in the 
direction that is followed by the beam inside the structure. This direction is given by the shift of the first transmitted 
beam (see figure 1), the shift being defined as the difference between the barycenter of the incident field and that of 

/ x\u° t (xfi)\ 2 dx 
f\u°(x,0)\' l! dx ' 

We assume that for the considered frequency all the plane waves inside the beam belong to the conduction band of 
the photonic crystal. This implies that I7I = 1, it can therefore be written under the form: 7 = e tky Nh . The vector 
(k x , k y ) is a "renormalized" Bloch vector in the following sense: denoting by u(x, 0) the value of the field on the lower 
interface, we have: u(x,Nh) = e lk:c x e lky Nh u(x, 0). 

Let us concentrate on the first transmitted beam, i.e. the beam that reads: 



ul (x, Nh) = J A(k x )(l - \g\ 2 )e ik * Nh e ik * x dk x (9) 
Using Parseval-Plancherel identity, we get the angular shift due to the beam propagation (cf. fig. 1): 

A I^(k x )(l-\g\ 2 ) 2 dk x 



Nh r ,0 „ s /, , ,2 N 2 



jA 2 (k x )(l 



(10) 
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A series expansion of A can be obtained provided the phase function is analytic with respect to k x in a neighborhood 
of k rn . Then we obtain: 



m=0 m! dkl 



m+l 



(11) 



wWp r - /^ 2 (M(iHg| 2 ) 2 (fc*-^o)'"^ 
wnere o m - j A 2 (fex)(1 _| g p)2 dfcx ■ 

When A (fc^) is concentrated around fc m , and if does not vary too quickly in the vicinity of k m , we obtain the 
well-know crude approximation: 

A ~ -Nh^L (k m ) (12) 

In order to give a geometric interpretation of this result, let us remark that {^f- {k m ) , —l\ is a vector that is 

normal to the dispersion curve at wavelength A. 

We retrieve the well-known fact that for a spatially large beam, the direction of propagation is given by the normal 



to the isofrequency Bloch diagram [15| but here it is the "renormalized" Bloch diagram that is involved. 

We will see in the numerical application that it can be quite different from the usual Bloch diagram. 

In the particular case of a one dimensional stratified medium, i.e. when the relative permittivity is constant in 
the horizontal direction, the direction of propagation of a beam is given by the normal to the usual dispersion curve. 
Indeed, denoting T the transfer matrix for 1 period, the transfer matrix for N periods is T N . This matrix coincides 
with the matrix T^r : this is due to the absence of evanescent waves. A direct consequence is that the Bloch vectors 
obtained from T N and T^v are the same, hence the renormalized Bloch diagram coincides with the non-renormalized 
one. 



DISCUSSION 



In the following, we present numerical computations illustrating the behavior of a beam inside finite PhCs. All 
the numerical results were obtained by means of a rigourous diffraction code for gratings based on the Fourier Modal 
Method (see [Hj]). We will denote: 

• A„ the shift of the first transmitted beam computed by a direct numerical computation of the field. 

• Ab the shift computed through the isofrequency Bloch diagram. 

• An the shift computed through the renormalized isofrequency diagram 

• Af the shift of the entire transmitted field comprising all the beams, computed numerically. 

In order to quantify the role of the evanescent waves, let us remark that, inside the photonic crystal, the field can be 
expanded over three types of modes [22| ■ These modes are the eigenvectors of the transfer matrix T of the photonic 
crystal, i.e. the matrix that relates the fields below the crystal to the fields above the crystal: 

1. the propagative modes, i.e. the Bloch modes (corresponding to the eigenvalues of T of modulus 1), 

2. the evanescent modes (corresponding to the eigenvalues of T of modulus less than 1), 

3. the anti-evanescent modes (corresponding to the eigenvalues of T of modulus greater than 1). 

By projecting the field on these modes, it is possible to compute the ratio of the field that is carried by the evanescent, 
anti-evanescent and propagating waves [22I ] . 

We consider a photonic crystal with square symmetry, infinite in the horizontal direction and comprising N periods 
in the vertical direction. The basic cell is given in fig. 2: it is made of square air holes (side d/2) inside a dielectric 
matrix of permittivity e = 9. We first compute the transmitted field for N — 1: the structure is made of one single 
grating. The incident beam is a gaussian p-polarized beam (w — 50d, X/d — 2.5, 9 = 40°). As it has been explained 
before the transmitted field is a sum of beams: the total field on the upper interface is given in fig. 3 and the 
amplitudes of the different beams \u^(x, Nh)\ for p = 0, 1,2,3,4 in fig. 4. It is clearly seen that there is a strong 
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overlap between the beams. The shift of the first beam is A n /d = 1.66 while the shift of the total transmitted beam is 
A f/d = 0.27. The Bloch diagram is given in fig. 5 (black curve). The value of k m is: k m xd = 2ir/2.5 x sin(40°) ~ 1.6. 
The predicted shift is Ag/d = —1.36, the renormalized Bloch diagram (grey curve in fig. 5) gives: A^/d = 1.66. In 
this situation, the overlap between the multiple beams is strong and the transmitted field cannot be reduced to the 
first beam. 

Let us now keep the same parameters except for the number of periods: we take N = 4. This time the situation 
changes drastically: the transmitted field shows clearly a negative refraction (fig. 6), we have A n /d ~ —36, the shift 
of the total field is A f/d ~ —32 and the Bloch diagram gives: As/d ~ — 11. In this situation the overlap between 
the beams is less important (fig. 7). The renormalized Bloch diagram gives: An/d ~ —27. The discrepancy between 
A n and A# is due to the fact that the beam is not spectrally narrow enough. If we use a narrower beam by taking 
w = 200c? we get a better result: A n /d ~ —28. Let us look at the modulus of the eigenvalues of T and the ratii 
on the different modes(fig. 8). On this figure, the modulus of the eigenvalues are given in thin solid line, the ratio 
over the propagative waves is shown in thick solid line and the ratio over the evanescent waves in dashed line. It is 
clearly seen that near X/d ~ 2.5 there are more evanescent waves than propagative ones: this explains why the Bloch 
diagram gives a false description. Let us pursue and take N = 16, keeping the same parameters as for TV = 4, we get 
a transmitted field given in fig. 9 the beam is hugely deformed and the value of its barycenter is not relevant: here 
the behavior of the field is much more complicated than could be expected from the Bloch diagram. 

Although the above results show that the prediction of the Bloch diagram can be quite false, it should be said that, 
provided the influence of the evanescent waves is not too strong, the Bloch diagram can give very accurate results. 
This is the case if we take X/d — 6 with the same parameters as above. We get: A n /d — 6.55 and Ag/d = 6.45. This 
time the Bloch diagram predicts fairly well the position of the first transmitted beam. 



CONCLUSION 



We have shown that the shift of the first transmitted beam is given by the normal to the renormalized isofrequency 
diagram and that it can be quite different from the prediction of the Bloch diagram. We have proven that this effect 
is due to the existence of evanescent waves inside the structure. The analysis uses a dressed transfer matrix that could 
be extended to aperiodic or random structures 23j, 124} . We have also pointed out the strong influence of the number 
of periods and of the spectral widths of the beams: these parameters can drastically modify the behavior of the beam. 
Moreover, it may happen that the total transmitted field be not properly described by the first transmitted beam. In 
that situation, neither the Bloch diagram nor the renormalized diagram can give reliable results. This study aims at 
showing that one should be prudent when using the Bloch diagram to predict the behavior of the field inside a true 
(i.e. finite) photonic crystal. Beyond that, it shows that many new effects can be expected due to the presence of 
evanescent waves. These are not a drawback but rather represent new possibilities to imagine new fonctionnalities. 
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Figures captions 

Figurel: 2D Photonic crystal as stacks of gratings: periodic in x and finite in y. 
Figure 2: Basic cell of the photonic crystal used in the numerical experiments. 
Figure 3: Amplitude of the transmitted field on the upper interface. 
Figure 4: Amplitude of the first five transmitted beams. 

Figure 5: Black line: Bloch diagram for the photonic crystal. Gray line: renormalized Bloch diagram. k x x d and 
k y x d belong to (— 7r, +tt). The vertical dashed line indicates the average value of k x in the incident field. 
Figure 6: Amplitude of the transmitted field on the upper interface. 
Figure 7: Amplitude of the three first transmitted beams. 

Figure 8: Thin solid line: modulus of the eigenvalues of the transfer matrix. Thick solid line (blue online): ratio of 
the field on the propagating modes. Dashed line (red online): ratio of the field on the evanescent modes. Dashed 
dotted line (green online): ratio of the field on the anti-evanescent modes. This arrows indicate the wavelengths were 
the computations were performed. 

Figure 9: Amplitude of the transmitted field on the upper interface. 
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